Revisiting implementation of multiple natural enemies in pest management

A major goal of biological control is the reduction and/or eradication of pests using various natural enemies, in particular, via deliberate infection of the target species by parasites. To enhance the biological control, a promising strategy seems to implement a multi-enemy assemblage rather than a single control agent. Although a large body of theoretical studies exists on co-infections in epidemiology and ecology, there is still a big gap in modelling outcomes of multi-enemy biological control. Here we theoretically investigate how the efficiency of biological control of a pest depends on the number of natural enemies used. We implement a combination of eco-epidemiological modelling and the Adaptive Dynamics game theory framework. We found that a progressive addition of parasite species increases the evolutionarily stable virulence of each parasite, and thus enhances the mortality of the target pest. However, using multiple enemies may have only a marginal effect on the success of biological control, or can even be counter-productive when the number of enemies is excessive. We found the possibility of evolutionary suicide, where one or several parasite species go extinct over the course of evolution. Finally, we demonstrate an interesting scenario of coexistence of multiple parasites at the edge of extinction.


SM1 Single Infection Model
The model equations for S and I when we only have one type of parasite present are as follows dS dt = F(N)S − µS − Sβ 1 1 I 1 , (S1) dI 1 dt = β 1 1 I 1 S − (µ + α 1 ) I 1 .
The stationary states of these are given as follows As we have F(N) = r(1 − N/K), This results in and N = α 1 + r β 1 1 + r/K .

SM1.1 Evolution with One Parasite Present
In this case our mutant equations are given as follows dI m dt = β m m I m S * − (µ + α m ) I m .
This defines invasion fitness as It is straight-forward to derive that S = µ+α 1 We can also use this expression to locate the ESS.
This gives the ESS for a single strain as α * s = Kµ. (S11)
A detailed flowchart of the model is provided in Figure S1. Figure S1. Flowchart displaying all interactions within the generic model that accounts for the existence of multiple strains of parasite, given by equations (S12)-(S14). The system consists of three components; healthy hosts (S), hosts infected by a parasite of strain i only (I i , i = 1, 2) and hosts infected by both strains of parasite (D).

SM2.1 Evolution with a Single Parasite
The model flowchart of this is shown in Figure S2, in terms of the adaptive dynamics framework, we consider S, I 1 , I 2 and D to be at the stationary states and therefore we just model the rare mutant strains of I 2 , denoted as I m and D 1m . As the mutant strains are so rare we consider I m D 1m ≈ 0 simplifying the invasion model to the following In matrix form this is given as follows which gives the following characteristic equation From which we can obtain λ 1,2 , we take the invasion fitness to be λ (α 2 , α m ) = max(λ 1,2 ). Figure S2. Flowchart displaying all interactions within the generic model that accounts for the existence of multiple strains of parasite, given by equations (S12)-(S14). The system consists of three components; healthy hosts (S), hosts infected by a parasite of strain i only (I i , i = 1, 2) and hosts infected by both strain of parasite (D).

SM2.2 Co-Evolution of Two Co-infecting Parasites
We also consider the case when we allow both strains of the parasite to evolve, i.e. co-evolution. The flowchart for this scenario is shown in Figure S3. In terms of the adaptive dynamics framework, we consider S, I 1 , I 2 and D to be at the stationary states and therefore we just model the mutant strains; Figure S3. Flowchart displaying all interactions within the generic model that accounts for the existence of multiple strains of parasite with the invasion of a mutant of both strains, given by equations (S19)-(S23). The system consists of three components; healthy hosts (S), hosts infected by a parasite of strain i only (I i , i = 1, 2, m 1 , m 2 ) and hosts infected by two strains of parasite (D i, j̸ =i ).
Again, note that due to how rare the mutants are I m 1 D m 1 = I m 1 D m 2 = I m 2 D m 2 = I m 2 D m 1 = 0. From this model we use λ m1 = max(λ 1 , λ 4 ) and λ m2 = max(λ 2 , λ 3 ) as the invasion fitness for I m 1 and I m 2 , respectively. As dD m 1 m 2 dt < 0 we can neglect it from out analysis.

SM3 Triple Infection Model
Now consider the case when we have three possible types of parasite; I i j (denoted here by I j , j = 1, 2, 3) which can pairwise doubly infect hosts; I i j ,i k (denoted here by D 12 , D 13 and D 23 ), furthermore, the host can be triply infected by all three types I i 1 ,i 2 ,i 3 (simply denoted by T ). The model for this scenario can be described by the following set of equations where

SM5 Cost-Effectiveness of Biological Control
In this study, we consider two measures of cost-effectiveness of biological control by co-infecting parasites: the average cost-effectiveness ratio (ACER) and the incremental cost-effectiveness ratio (ICER), for details of the formal definitions, see the main text. Figures 3 (C,D) of the main text show the values of ACER. In figure S4, we show the graphs of the ICER corresponding to the parameters in Figure 3 of the main text. The negative value for the ICER observed for n = 4 in panel (A) signifies a decrease in the effectiveness with an increase in the cost while transferring from the triple to the quadruple co-infection biological control. Finally, we analyse the cost-effectiveness of biological control according to the scenario, where the efficient population size corresponds to the total number (density) of the susceptible host, i.e. N e = S. In this case, the effectiveness is measured as (S(0) − S(n))/S(0), and the increment in effectiveness -when adding an extra type of natural enemy in the system to the existing n − 1 types -is defined as (S(n − 1) − S(n))/S(0). The corresponding graphs for both the ICER and the ACER are shown in Figure S5. Figure S4. Dependence of the ICER on the possible number of co-infections n for two cases of transmission β i j i 1 ,i 2 ,...,i k . The blue curves display the ICER when effectiveness is evaluated using the total species density N at the co-ESS virulence (mathematically, the increment in effectiveness, when adding an extra type of natural enemy to the existing n − 1 types is given by (N(n − 1) − N(n))/N(0)). The orange and yellow curves display the ICER based on the efficient population size of the pest N e at the ESS virulence (mathematically, the increment in effectiveness, when adding an extra type of natural enemy to the existing n − 1 types, is given by (N e (n − 1) − N e (n))/N e (0)), when weighting functions w 1 and w 2 are used. The dashed curves correspond to the short-time scale scenario, with the virulence given by the ESS for a single infection (n = 1). Within the weighting functions w 1 and w 2 , the parameters are as follows; α 0 = 0.4, a = 1, b = 20 and c = 5, all other model parameters are as defined in Table 1. Figure S5. Dependence of (A) the ICER and (B) the ACER on the possible number of co-infections n for two scenarios of transmission β i j i 1 ,i 2 ,...,i k in the case when the effectiveness is evaluated based on the total susceptible host density S at the co-ESS virulence. The dashed curves correspond to the short-time scale scenario, with the virulence given by the ESS for a single infection (n = 1). Within the weighting functions w 1 and w 2 , the parameters are as follows; α 0 = 0.4, a = 1, b = 20 and c = 5, all other model parameters are as defined in Table 1.

SM6 Calculating the prevalence of co-infections and the corresponding mortality rates.
In this section we evaluate the values co-infection prevalence p(I i 1 ,i 2 ,...,i k ) (i j = 1, ..., n) for the system at equilibrium. To find the equilibrium values of host density, we use the model simulation in the main text. For each maximal number n of possible co-infections, the prevalence of a host having a particular number of co-infections is determined by the following proportion where N is the total number of hosts, in the above expression below we understand all densities as stationary ones. Table 1 presents the co-infection prevalence values corresponding to the scenarios shown in Figure 6 of the main text, where S, I, D, T, Q denote, respectively, the total population size (the total density) of the host infected by k = 0, k = 1, k = 2, k = 3, and k = 4 different types of parasites with k ≤ n, n being the total maximal number of co-infections (e.g. I = I 1 + I 2 + ...I n ).  Table 1. Prevalence (proportion) of host containing various numbers of co-infections (n the maximal number of co-infections), for the two different cases of transmission rates, calculated based upon the stationary densities displayed in Figure 6 of the main text.

9/10
In this section, we also calculate the probability that a host will die after having acquired exactly n infections (for generality, we consider n = 0, i.e. no infections, as well). For the system at equilibrium, mentioned probabilities can be calculated based on the following expression p m (I i 1 ,i 2 ,...,i k ) = (µ + α i 1 ,i 2 ,...,i k )I i 1 ,i 2 ,...,i k µS + ∑(µ + α i 1 ,i 2 ,...,i k )I i 1 ,i 2 ,...,i k = (µ + α i 1 ,i 2 ,...,i k )I i 1 ,i 2 ,...,i k F(N)S , where µ is the background mortality, α i k 1 ,i k 2 ,...,i kn is the virulence due to infection by parasites i k 1 , i k 2 , ..., i k n . The rationale behind the above expression is that for the system at equilibrium, the reproduction rate (inflow) of the host should be compensated by the losses due to mortality. On the other hand, the number of hosts dying after acquiring particular infections i k 1 , i k 2 , ..., i k n is proportional to the magnitude of the corresponding mortality rate, i.e. it is proportional to (µ + α i 1 ,i 2 ,...,i k )I i 1 ,i 2 ,...,i k . Importantly, the above expression defines the probability that a host will be infected by the maximal number of co-infections n, i.e. dies after acquiring all possible types of co-infections: p m (I 1,2,3,...,n ) = (µ + α 1,2,3,...,n )I 1,2,3,...,n F(N)S .
Using the above equations for p m , we calculated the probabilities to die for host individuals as a result of acquiring exactly k infections, k ≤ n. The results are presented in Table 2 which corresponds to modelling graphs shown in Figure 6 of the main text.  Table 2. Probabilities of host depth after having acquired k infections, k ≤ n (k = 0, 1, 2, 3, 4 correspond to S, I, D, T, Q, respectively). Two different cases of transmission rate are considered, the stationary densities of species are the same as displayed in Figure 6 of the main text.
A comparison of Table 1 and Table 2 shows that the probability to die after k co-infections p m can be substantially different from the corresponding values of co-infection prevalence p for the same number of co-infections. This is especially true for the maximal number of co-infections k = n. The mentioned discrepancy can be explained by an increasing mortality rate of hosts with a large number of co-infections.